% function to calculate bivariate errors
% given coherency and spectra
% coh - e-pred coherency for current eqn
% df degree of freedom
% data spectral matrix
% Ref: ProcMT manual
% not yet tested OK
% tested OK ! 02.08.2002
% falg 'Ex' or 'Ey' (indicating eqn)

% latest date 3/5/2003


function[da,db] = BivError(coh, df, data, flag),

j = sqrt(-1);
coh=coh^2; % squared predicted coherence added 3.5.3
HxHx = data(1,1);
HyHy = data(2,2);
HzHz = data(3,3);
ExEx = data(4,4);
EyEy = data(5,5);


HxHy = data(2,1) - j*data(1,2);
HyHx = conj(HxHy);
HxEx = data(4,1) - j*data(1,4);
ExHx = conj(HxEx);
HxEy = data(5,1) - j*data(1,5);
EyHx = conj(HxEy);
ExHy = data(4,2) + j*data(2,4);
HyEx = conj(ExHy);
HyEy = data(5,2) - j*data(2,5);
EyHy = conj(HyEy);



if flag == 'Ex',
ZZ = ExEx;
elseif flag == 'Ey',
ZZ = EyEy;
else,
disp('Only Ex or Ey allowed. Error ');
return;
end;



fis = fischer(df-4); % calculate fischer distribution F (4,df-4,0.05) Bendat & Piersol, page 110

d =  (4/(df-4)) * fis * (1.0 - coh) * ZZ;

if abs(HxHx)*abs(HyHy) > 0,
r2xy = (HxHy*HyHx)/(HxHx*HyHy);
else,
r2xy = 100;
end;

if r2xy > 0.999,
   r2xy = 0.999;
end;

if abs(HxHx*(1-r2xy)) > 0,
da = d / (HxHx*(1-r2xy));
else,
da = 100;
end;

if abs(HyHy*(1-r2xy)) > 0,
db = d / (HyHy*(1-r2xy));
else,
db = 100;
end;

da = sqrt(real(da)); 
db = sqrt(real(db));